#ifndef _SYSTEM_HF_H_
#define _SYSTEM_HF_H_
#include "hfsolver_ho.h"
#include <eigen3/Eigen/Eigen> //matrix library
#include <stdio.h>

/*
a common interface for method based on but beyond HF method.
System_TzType and DataType can be changed according to needs.
*/

using Eigen::Dynamic;
using Eigen::Matrix;
using Eigen::SelfAdjointEigenSolver;

typedef double HFDataType;
/// system under HF basis
class System_Tz_HF : public System_Tz<SHF_Orbital<HFDataType>, HFDataType> {
public:
  typedef System_Tz_HO System_TzType;
  typedef SHFSolver_HO SHFSolverType;
  typedef typename System_TzType::DataType DataType;
  /* typedef System_Tz_WS System_TzType; */
  /* typedef SHFSolver_WS SHFSolverType; */
  /* typedef typename System_TzType::DataType DataType; */
  typedef Matrix<DataType, Dynamic, Dynamic> Mat;
  typedef Matrix<DataType, Dynamic, 1> Vec;

  typedef SHF_Orbital<DataType> OrbitalType;
  const static int MatEledim; //< dimension of TwoBodyMatEle.
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  typedef Matrix<TwoBodyMatEle, Dynamic, Dynamic> Mat2B;
  int Flag_TBME = 0; // ! changed by zhuo
                     // 0: get2BmatOnNeed
                     // 1: setupTwoBodyMat()&get2BmatAt
                     // 2: get2BmatAt only, For NNlosat

  System_Tz_HF() : beta(0) {}
  System_Tz_HF(const System_TzType *Ptr)
      : pSystem_Tz(Ptr), hfsolver(Ptr), hbar_omega(Ptr->hbar_omega), beta(0) {}
  System_Tz_HF(const System_TzType *Ptr, int _A, int _Z)
      : pSystem_Tz(Ptr), hfsolver(Ptr, _A, _Z), A(_A), Z(_Z), realA(_A),
        hbar_omega(Ptr->hbar_omega), beta(0) {}
  System_Tz_HF(const System_TzType *Ptr, int _A, int _Z, int Flag_TBME_t)
      : pSystem_Tz(Ptr), hfsolver(Ptr, _A, _Z), A(_A), Z(_Z), realA(_A),
        hbar_omega(Ptr->hbar_omega), Flag_TBME(Flag_TBME_t), beta(0) {}

  // read interaction under HF basis directly from file
  void setupFromFile(const char *orbfile, const char *intfile);

  void setAZ(int _A, int _Z, int _realA = -1) {
    if (_realA == -1)
      realA = _A;
    else
      realA = _realA;
    A = _A;
    Z = _Z;
    hfsolver.setAZ(_A, _Z, realA);
  }

  void setupOrbitals();
  void setupTwoBodyMat(); 

  void setup(const System_TzType *Ptr, int _A, int _Z) {
    pSystem_Tz = Ptr;
    hbar_omega = Ptr->hbar_omega;
    beta = 0;
    A = _A;
    Z = _Z;
    realA = _A;
    hfsolver.setup(Ptr, _A, _Z);
    setup();
  }
  void setup() {
    setupOrbitals();
    setupTwoBodyStates();
    if (Flag_TBME == 1)
      setupTwoBodyMat();
  }
  /// interface for method beyond HF.
  DataType get1B(int bra, int ket) const;

  DataType get2B(int a, int b, int c, int d, int J) const {
    TwoBodyMatEle MatEle;
    if (Flag_TBME == 1 || Flag_TBME == 2) {
      get2BmatAt(a, b, c, d, J, MatEle);
      // !! zhuo get directly also need setupTwoBodyMat()
      // !! Or nnlosat
    } else if (Flag_TBME == 0) {
      get2BmatOnNeed(a, b, c, d, J, MatEle); //
    } else {
      cout << "Wrong at system_hf.h :  get2B" << endl;
      exit(0);
    }
    if (MatEle.empty())
      return 0;
    if (a == b)
      MatEle *= sqrt(2.0);
    if (c == d)
      MatEle *= sqrt(2.0);
    return MatEle[0] - MatEle[3] * pSystem_Tz->hbar_omega / realA;
  }

  /// interface for nocore shell model
  DataType get1B_p2(int bra, int ket) const;
  double get1B_p2_plus_r2(int bra, int ket) const;

  void get2BmatOnNeed(int a, int b, int c, int d, int J,
                      TwoBodyMatEle &MatEle) const;

  void printOrbitals(int num = -1) const {
    if (num == -1)
      num = Orbitals.size();
    cout << "num"
         << "\t"
         << "l"
         << "\t"
         << "jj"
         << "\t"
         << "tz"
         << "\t"
         << "e"
         << "\t"
         << "OccuP" << endl;
    for (int i = 0; i < num; i++) {
      cout << i << "\t" << Orbitals[i].l << "\t" << Orbitals[i].jj << "\t"
           << Orbitals[i].tz << "\t" << Orbitals[i].e << "\t"
           << Orbitals[i].OccuP << endl;
    }
  }

  void setbeta(double _beta) { beta = _beta; }

  /// maybe usefull operator
  /// only radial part
  DataType rL(int bra, int ket, int L) const;
  DataType rms() const {
    DataType val = 0;
    for (int i = 0; i <= FermiSurface; i++) {
      val += rL(i, i, 2) * (Orbitals[i].jj + 1.);
    }
    return sqrt(val / realA);
  }

  DataType rms_p() const {
    DataType val = 0;
    for (int i = 0; i <= FermiSurface; i++) {
      if (Orbitals[i].tz == -1)
        val += rL(i, i, 2) * (Orbitals[i].jj + 1.);
    }
    return sqrt(val / Z);
  }

  double HF_r(int orb, double r) const {
    double temp = 0;
    for (int i = 0; i < Orbitals[orb].Coeff.size(); i++) {
      int iorb = pSystem_Tz->Groups[Orbitals[orb].GroupIndex][i];
      temp += pSystem_Tz->HO_r(iorb, r) * Orbitals[orb].Coeff[i];
    }
    return temp;
  }
  /*
    Q0=sum_i{r_i^2},S(E0)=2\hbar^2/m*A*<r^2>
  */
  /// be careful with the convention, Q0 is r^2, while QL proportion to r^L YL
  DataType Q(int p, int h, int J, int isospin) const;
  DataType Q1_v(int p, int h) const;
  DataType Q1_s(int p, int h, DataType Rms) const;

  /// type=1 for T_plus, -1 for T_minus
  //<a| Beta_F |b> = <a| t+ or t-| b>
  DataType Beta_F(int type, int a, int b) const;
  //<a|Beta_GT|b> = <a|| sigma * (t+ or t-)  ||b>
  DataType Beta_GT(int type, int a, int b) const;
  DataType radialOverlap(int a, int b) const;

  int A, Z;  // for HF. basis generating.
  int realA; // the real A in Hamiltonian which remove COM.
  int FermiSurface;
  DataType HFenergy;
  const System_TzType *pSystem_Tz;
  SHFSolverType hfsolver;

  double beta;
  double hbar_omega;

private:
  /// return <a|alpha>, a is the ath orbital of Orbitals in pSystem_Tz, alpha is
  /// alpha'th HF. orbital in Orbitals
  /* DataType overlap(int a,int alpha) */
  /* { */
  /*   if(Orbitals[alpha].GroupIndex != pSystem_Tz->GroupIndice[a]) */
  /*     return 0; */
  /*   return Orbitals[alpha].Coeff[pSystem_Tz->OrbIndice[a]]; */
  /* } */
};

// for norm order nnlo_sat interaction
void setup_nnlo_sat(System_Tz_HO &system_ho, System_Tz_HF &system_hf,
                    const char *spfile, const char *umatfile,
                    const char *kinfile, const char *vfile);
double nnlo_sat_getHFenergy(const System_Tz_HF &system_hf);
#endif
